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SUMMARY 


We show that by certain transformations the boundary layer equations for the class of 
non-Newtonian fluids named pseudoplastic can be generalized in the form 

A u + p{x)u~ x = 0, ieOC R n , n > 1 

under the classical conditions for steady flow over a semi-infinite flat plate. We provide a survey of 
the existence, uniqueness, and analyticity of the solutions for this problem. We also establish 
numerical solutions in one- and two-dimensional regions using multigrid methods. 

INTRODUCTION 

In the last two decades, solutions of the singular semilinear equation 

A u + p(x)u~ x — 0, x E Q C R n (1) 

have been extensively studied. Various existence and uniqueness results are given in [1], [2], and [3], 
to name a few. More recently, in [4], it is shown that by certain transformations the boundary layer 
equations for the class of non-Newtonian fluids named pseudoplastic can be generalized in the above 
form for the ODE case n = 1. Under this physical interpretation the above equation, considered in 
the context of partial differential equations (n > 1), has been the subject of much study. The 
equation has a unique classical solution with a bounded domain fi, where p{x) is a sufficiently 
regular function which is positive on Q. [5]. There exist entire solutions with A G (0, 1) for p(x) 
sufficiently regular ([6], [7]). This is generalized to all A > 0 via the upper and lower solution 
method ([8]) or other methods ([9]). 
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The following sections provide a survey of both theoretical and numerical results in this area 
including a physical derivation [4], existence theorems for both the ODE and PDE cases with a 
proof of a main result [8], and our numerical results. We conclude with a discussion of a new 
technique and some open questions for further research. 

PRELIMINARIES 


A non-Newtonian fluid is called pseudoplastic if the shear stress r and the strain rate 



related as 


|r| = k 


du 

dy 


Of 


J 


0 < a < 1 


where k is a positive constant. That is, the absolute value of the shear stress increases with respect 
to the absolute value of the strain rate less than linearly. 

In this paper, we study solutions of the singular semilinear equation (1) where A > 0 and Q is a 
domain in R n , n > 1. In the following section we show that through a series of transformations the 
boundary layer equations for the class of pseudoplastic fluids under the classical conditions for a 
steady flow over a semi-infinite flat plate can be generalized into the well-known Blasius problem 


+ //' = 0 , 


f{0) = /'(0) = 0, /'( oo) = l 


for the shear function, which arises from the standard Newtonian fluid case. 


DERIVATION OF THE PROBLEM 


For n = 1 equation (1) arises in the study of pseudoplastic fluids. We consider a two-dimensional 
incompressible flow of low viscosity along a plane wall. We denote by v = ( u , v) as the fluid velocity 

in the boundary layer and u. x (x) the main stream. Since there is no velocity on the wall and the 

du 

fluid takes the velocity of the main stream u x .(x) outside the boundary layer, we see that is 

large near the wall which causes a significant transfer of momentum in the x direction. 

The boundary layer equations for this model include a continuity equation and a momentum 
equation in the x direction. 


du du _ 
dx + dy 


u 


du du 1 dr. 


dx 


+ v— = -• 
dy p 


xy 


dy 


( 2 ) 


with boundary condtions 
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it(x,0) = v(x,0) = 0, 
u(x,oo) = u. x (x) 
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Figure 1: 2-D flow of low viscosity along a plane wall. 


where r iy = K 


du 


dy 


is the shear stress. 


Note that (2) has 2 coupled equations in 3 dependent variables, u, v, and r xy . To reduce it to a 
single higher-order equation in only 1 dependent variable, we introduce the Lagrange stream 
function <h(x, y) such that 

Then the momentum equation becomes 

a$ a 2 # a$ a 2 $ 


dy dxdy dx dy 2 dy K dy 2 


( 3 ) 


K 


where v = — , while the continuity equation is clearly satisfied by <I>. 

p 


Let / = 


y/uf~ 


i/x 


; r) = byJ — for some a, b. Then (3) becomes 
v vx 

r+/(n 2 -“=o 


( 4 ) 


with 


/( 0 ) = /'( 0 ) = 0 

/'(°°) = 1 

where / = f(f]). (Observe that if a = 1 (4) is the well-known Blasius equation.) Employing the 
Crocco-like transformation 

^ = /'(*?), g(u) =v! = f'fn) 


(4) becomes 


g a g" + (a-l)g a - 1 (g , ) 2 + u = 0 


with y'(0) = 0, 5(1) = 0, where g = g(u). Finally the transformation G = g a leads to the singular 
boundary value problem 

G" + auG~ l/a = 0, 0< u,a<l, 

G'{ 0) = G(l) = 0 

1 X 

of the form (el) with A = — , u — x, P = ~r- 


EXISTENCE AND UNIQUENESS RESULTS 


In the first part of this section we study the results in finite and infinite domains; in the second 
we discuss methods that are commonly used to approach the problem. 


Theorems 


Let Fl be a bounded domain in R n , n > 1 with smooth boundary dQ (of class C 2+Q , 0 < a < 1). 
Let p(x ) be of G“(D) and positive on Cl, A > 0. 

Theorem 1 (Lazer-McKenna [5]). The problem 

A u + p(x)u~ x =0, 
u |sn= 0 

has a unique positive solution u(x) in Q with u 6 C’ 2+ “(D) D (D). Furthermore let be an 
eigenfunction corresponding to the smallest eigenvalue Ai of the problem 

A<f) + \(p = 0, x eft 
4> |an= 0 

such that ef > i(x) > 0, 1 6 1 ) and A > 1. Then there exists a unique 61,62 > 0 such that 

M; /(,+A) < U < W? /<1+A) 

on D. 


In the case Q = R n , n > 1, we study the results under conditions n = 1, n = 2, n > 3. Observe 
that if n = 1, since p, y > 0, y” + py~ x = 0 we have y" > 0 and thus y' j. Hence 0 < y'( 00 ) < 00 . 

Theorem 2 (Taliaferro [3]) The problem 

y" + p(x)u~ x — 0 
y(c) = a 

y '{ 00) = 0 

has a unique positive solution y(x) if 

j x~ x p[x) dx < 00 

where a, c £ R 1 , a > 0. Furthe rmore y( 00 ) <00 if and only if ff' xp(x)dx < 00 . 
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The following theorem describes the asymptotic behaviors of the solution. 

Theorem 3 (Taliaferro [3]) 

• If 0 < y'( oo) < oo and Jj 00 x~ x+1 p(x)dx < oo, a, b > 0 then 

y{x) — ax + b — a-A(l + o(l)) f (f - x)f _A p(f)d£. 

Jx 

• If y’( oo) = 0 and Jq° xp(x)dx < oo. a > 0 /hen 

y(x) = a - a— A(1 + o(l)) f (£ - x)p(£)e££- 

Jx 

• if p,q > 0 air continuous on [0, oo). lim x _ = R > 0 and 

2 " + p(x)z~ x = 0, z'(oo) — 0; 
w" + q(x)w~ x = 0, w'(oo) = 0 

and /» xp{x)dx — 00 , i/ten lim I _, x 10/2 = /fT* _ 

Theorem 4 (Kusano-Swanson [7]). The problem 

Au = /(|x|)it _A = 0, x € .R 2 , 0 < A < 1 

has an entire positive solution in R 2 with logarithmic growth at 00 if f(t) > 0. t > 0, 
f(t) £ C(0, 00 ), and 

J t(\ogt)~ x f(t)dt < 00 . 

A function u(x) is said to be an entire solution of (1) if u E Cf oc (R n ) and u satisfies the equation 
pointwise in R n . 

Theorem 5 (Shaker [8]) The problem 

Ait p(x)u~ x =0, x £ R n , A > 0 
has an entire positive solution it(x) such that 

ci < it(x)|x| 9|n - 21 < c 2 
for some C; , c 2 and 0 < q < 1 as x — ► OO if 

1. p{x) E Ct(R n ), P(x) > 0 for x E R n \{ 0> 

2. the re exists 0 < c < 1 such that cej)( |x|) < p(x) < 0(|x|) where ef>(t) = max| x | = t p(x) . t E [0, oo), - 

3. f? t n - 1+ ^ n ~ 2 ^(t)dt < 00 . 
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Methods 


In general there are two methods that are commonly used in proving existence and uniqueness of 
solutions for equations of type (1), namely Schauder’s fixed point thocrem and Barrier Methods. 
Since the former is standard we elaborate here only on the latter. 

Let Q be a smoothly bounded domain in R n . fi(x) is said to be an upper (lower) solution of the 
problem 

An + f(x, u) = 0, x G Q (5) 

w |en= 0 

if A <f> + /(x, fi) < 0, x G f2, fi(x) >0x6 dLl (fi + /(x, fi) > 0, x € fl, <j>(x) < 0 x G dSV). 

Theorem 6 (Sattinger [10]). Let <f>\ be an upper solution and (f> 2 be a lower solution oj (5), and let 
f be locally Holder continuous in Q. If fi i(x) > fii(x) in Q, then (5) has a solution u such that 
(f> 2 (x) < u(x) < fix (x). x G LI. 

In the case when Q. — R 71 we say fi is an upper (lower) solution of 

Au + f(x,u) = 0 (6) 

if A <f> + /(x, fi) < 0 x G R n (for lower solution, Afi + /(x, fi) >0). 

Theorem 7 (Ni [11])- Let fix and <f >2 be an upper and a lower solution of expiation (6). such that 
fix > <h x G R n ■ If f is locally Holder continuous in x and locally Lipschitz continuous in u. then 
(6) has a solution u with fii ( 2 ) < u(x) < fix (x) . x G R n ■ 

An Example. 

Consider the problem 

u" + Xu — it 3 = 0, x G (0,7r) 

U — 0, X = 0, 7T. 

It is easy to show that fix(x) = Rx l/2 for some R large is an upper solution, and fi 2 (x) = esmx for 
some e small is a lower solution of this problem. Clearly fi\(x) > fi 2 (x) for x G [0,7r]. Thus by the 
above theorem there is a solution u(x) such that e sin x < u{x) < Rx^ , x G [0 , xr] . Since the problem 
is homogeneous we conclude that the problem has at least three solutions, namely, u, —u and the 
trivial solution. 


MULTIGRID SOLUTION OF THE PROBLEM 


In this section we present some numerical results for solving the problem 

Ait + p(x) u~ x = 0 x G n 
u(x) = 0 x G dLl. 
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Specifically, we describe Newton’s method for non-linear systems to solutions and multigrid 
V — cycle and FMV methods. We have implemented all of these methods for both the one- and 
two-dimensional cases, using (respectively) the unit interval and the unit square for fi. In each case 
we use a straightforward finite-difference discretization, employing the standard second-order 
difference approximation for the second derivative operator. For the one-dimensional problem we 
define the grid of ( N + 1) points x k = jh, for k = 0, 1, . . . N, where h is the mesh parameter l/N. 
The second derivative operator is then approximated by 

= + . “*+■ + 0 (h% (7) 


dru 

dx 2 


where it* approximates ti(xjfc). For the non-linear term p(x)u(x) x we use the nodal values, pkU k x . 
Since u 0 = Un = 0, this results in the non- li near system of equations 


'-2 1 


u l 


p 1 u 1 x 


' 0 ' 

1 -2 1 


1*2 


P 2U2 X 


0 




+ 

: 

= 

0 

1 -2 1 


Un-2 


PN- 2Ujf X _ 2 



1 - 2 . 


l UN - 1 


p N -\ Uji * j 


0 


( 8 ) 


Letting u represent the vector of unknowns, we may write the system as Hu + g{u) = 0, where H 
is the tridiagonal matrix and g is the non-linear vector function. 

For the two-dimensional case we take the tensor product of the ( N + l)-point grid in the x 
direction with an identical ( N + l)-point grid in the y direction, yielding an ( N 4- l) 2 -point regular 
grid covering the unit square. The difference operator for the two-dimensional problem is 


( 


d 2 u d 2 u' 

d^ + d^ / 


' x j,k 


rr 


Numbering the unknowns lexicographically by lines of constant y 7 we obtain the nonlinear system 


' A B 


u i 


W\ 


' o ' 

BAB 


**2 


w 2 


0 

• . 


| 

+ 

; 

= 

0 

BAB 


UN- 2 


WN- 2 


j 

B A 


Ujv-1 


W N - 1 


0 


( 10 ) 


where here Uj denotes the ( N — l)-length vector of unknowns Uj k for k = 1, 2, . . . , JV — 1 
corresponding to the j th grid-line in the y direction, and A and B are (N — 1) x (N — 1) matrices 


'-4 1 


-1 

1 -4 1 

1 

-1 

* - 

B = — 

• . 

1 -4 1 

h 2 

-1 

1 -4 


-1 
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The (N - l)-length vectors Wj contain the non-linear entries p j}k Uj X k , for k = 1, 2 , . . . , JV - 1. Once 
again, we may write the system as Hu + g{ u ) = 0, where H is the block tri-diagonal matrix and g 
is the non-linear vector function containing the Wj's. 


Solution techniques 


The classical solution technique for (8) or (10) is to apply Newton’s method for non-linear 
systems. We write the system as F(u ) = 0, where F{u ) = Hu + g{u). Each step of the iteration 
is then given by 

u < — u — [J f ('u)] _1 F(u) 
where the Jacobian of the system is given by 

= H + D 

with H the linear part of F and D a diagonal matrix whose diagonal entries are the derivatives of 

the entries of g, for example -A . 

Naturally, the Jacobian is not inverted at each step, but rather, we solve the system 
[J/r(u)]y = —F(u) and then make the correction u <— u + y. We examined two methods for 
solving the system at each step, namely LU decomposition and a multigrid FMV cycle. 

Newton’s method converges quadratically. However, since each step involves inverting a system, 
it tends to be very slow. While the use of the FMV solver speeds the method up somewhat, it still 
is slower than the techniques we present next. It has long been known ([12], [13]) that on certain 
problems non-linear analogs to the classical Jacobi or Gauss-Seidel iteration methods could be 
employed with some success. Technically, one sweep of such a method means that for 
j = 2, . . . , N — 1 (or (N — l) 2 for the two-dimensional problem) one solves, via the scalar 

Newton’s method, the j th non-linear equation in the system F(u) = 0 for the j th unknown. As in 
the linear case, if the old values u are used throughout the sweep this is the Newton-Jacobi method, 
while if the updated values are employed as they become available it is the Newton-Gauss-Seidel 
method. In practice the j th equation is not actually solved, but rather, a few (one or two) steps of 
the scalar Newton’s method is performed on each equation in turn. 

The Newton-Jacobi and Newton-Gauss-Seidel iterations, however, typically behave in the same 
fashion that is observed in their linear counterparts. That is, the iteration generally progresses 
rapidly toward a solution with the first few sweeps, but then stalls out so that each additional sweep 
produces very little improvement. The reason behind this is the same as that seen in the linear case. 
The method stalls after the non-linear relaxation has successfully eliminated the oscillatory portion 
of the error, which it eliminates rapidly, but is unable to effectively treat the smooth portion of the 
error. This is precisely the difficulty that multigrid methods were devised to overcome. 

At the heart of multigrid is the coarse-grid correction [14]. Many common relaxation iterative 
relaxation methods for solving a linear problem Au = f have the property that the relaxation 
effectively eliminates the high-frequency (oscillatory) components of the error but leave the low 
frequency (smooth) components essentially unaffected. However, because the error is smooth after 
the relaxation, it may be represented accurately on a coarser grid, on which it also appears more 
oscillatory (relatively). Relaxation on this coarser grid then eliminates the oscillatory components 
of the coarse-grid error, which cannot be eliminated on the fine grid. The coarse-grid correction for 
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( 11 ) 


a linear problem may be written as 

u h «- P v u h + l 2 h (A 2h y 1 ll h (f h - A h P v u h ) 

where P is the relaxation matrix, v is the number of relaxations, is a prolongation or 
niU rpolahon matrix mapping coarse-grid vectors to the fine grid, l\ h is a restriction matrix 
mapping fine-grid vectors to the coarse grid, and A 2h is a coarse-grid version of the original matrix 
A. A crucial feature is that on the coarse grid Q 2h , the problem to be solved is the residual equation 
At = r, where the residual is defined r — f — Au and e is the error. That is, if u* is the exact 
solution, then At — A{u* — u) — f — Au = r. 

For nonlinear problems the residual equation doesn’t hold. Instead, we write the nonlinear 
equivalent of the residual equation, 

F(u + e) — F(u) = t. 

This equation is to be solved on the coarse grid, so we write 

F 2h {I 2 h h u h + t 2h ) - F 2h {I 2 h h u h ) = I 2h (f h - F h (u h )), (12) 

or 

F 2h (u 2h ) = I 2h (f h - F h (u h )) -I- F 2h (I 2h ). 

The coarse-grid correction is then performed by solving (12) for u 2h = ll h u h + t 2h , and then 
making the correction u h *— u h + /^('U 2 ^ — I^u h ). This gives the full approximation scheme [15] 

u h *- P> A ) + I% h {{F 2h )-\l 2 h h {f h - F\P v {u h ))) + F 2h {lfP v {u h ))) - I 2 h h P u (u h )), 
where P is a nonlinear relaxation scheme. 

For both the linear and nonlinear problems, the solution of the coarse-grid problem is computed 
using the same coarse-grid correction scheme as is being employed to solve the fine-grid problem. 
This leads to the multigrid V-cycle scheme, which (for the nonlinear problem using FAS) is 
described recursively as follows. 

«*«- FASV h (u h , f h , i/j, i/ 2 ) 

1. Perform Ui non-linear relaxation sweeps times on F h (u h ) = f h with initial guess u h . 

2. If Q h is the coarsest grid, then go to 4. Else: 

f 2 h _ jihyh _ F h( u h)) + F 2h (I 2h u h ) 
u 2h *— 0 

u 2h FASV 2h {u 2h ,f 2h ,v 1 ,v 2 ). 

3. Correct u h u h + I% h (u 2h — ll h u h ). 

4. Perform u 2 non-linear relaxation sweeps times on F h (u h ) = f h with initial guess u h . 

An important consideration for this (or any) iterative method is the choice of a good initial 
guess. Clearly a better initial guess will reduce the overall effort required to obtain an acceptable 
solution. A standard approach in multigrid is to obtain a good initial guess by first solving the 
problem on a coarse grid, and then interpolating that solution to the fine-grid for use as an initial 
guess. Solving this coarse-grid problem, in turn, will be easier if an initial guess is obtained by first 
solving the problem on a still coarser grid. Applying this idea recursively leads the Full Multigrid 
FMG scheme, which (applied to the non-linear / I S 1 scheme) may be described as follows: 

u h <- FASFMG h (u h , i /\ , v 2 ) 
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1. If is the coarsest grid, then go to 3. Else: 

f 2h = j2h(fh _ F h( u h )) + F ^^j2h u h^ 

u 2h «— 0 

u 2h <- FASFMG 2h {u 2h , f 2h , v x , U 2 ). 

2. Correct u h <— u h + I 2 hU 2h - 

3. u h *- FASV h {u h ,f h ,v x ,v 2 ). 


Numerical results for multigrid methods 


We have implemented the F AS V using Newton- Jacobi and Red-black Newton- Gauss- Seidel 
iteration schemes. (Our implementation was in Matlab using vector arithmetic. We elected not to 
analyse Newton-Gauss-Seidel since it is not vectorizable. We did encode it, however, and found 
that the slowness of the for loops overwhelmed the speed of convergence.) Several different choices 
for A, p(x) and p(x, y ) were used, as were several sets of relaxation parameters. 

Table 1 gives some quantitative information regarding the performance of the method, 
comparing convergence rates for various choices of parameters. The results shown were obtained 
using the Red-black Newton-Gauss-Seidel relaxation. We find that for this problem we are able to 
obtain convergence rates that are similar to those obtained on the linear elliptic model problems for 
which multigrid is best known ([14], [16], [17]). Data for the one-dimensional problem are not 
shown, however, they are very similar to the two-dimensional case. 


Dimension 

p{x) 

A 

Fine-grid 

Average V-cycle 



size 

convergence factor 

2 

2xy 

2 

63 x63 

0.051 H 



5 


0.050 



8 


0.078 

2 

2sin(27rx) sin(7ry) 

2 

63 x 63 

0.060 


5 


0.063 



8 


0.104 

2 

xjy 

2 

63 x 63 

0.059 


2 


0.060 



8 


0.086 


Table 1 


Additionally, we have implemented the F ASF MG using Newton- Jacobi and Red-black 
Newton-Gauss-Seidel iteration schemes. Again, we find that the performance of the method is 
compatible with that found for FMG applied to the linear model problems ( [15], [17]). 

CONCLUSIONS 


Our survey of existence and uniqueness results has shown the problem 

A u + p(x ) u~ x =0 x G Q 
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is guaranteed to have unique solutions under certain conditions, although these solutions will not be 
known in closed form. The problem arises in certain non-Newtonian fluids problems, so there is 
some interest in actually computing solutions. We have shown that for homogeneous Dirichlet 
boundary conditions on the unit interval and the unit square, multigrid methods appear to provide 
an efficient means of solution for reasonable choices of p(x). 

We note, however, that an actual convergence proof for the FAS method would be very difficult 
to obtain, in that such proofs normally require that we be able to decompose the space of grid 
functions into two operator-subspaces. Error components in one are annihilated by relaxation, 
while those in the other subspace are annihilated by coarse-grid correction. While such analysis is 
achieved for linear problems, non-linear problems generally can only be treated by linearization 
near a solution. In point of fact, the literature is remarkably sparse in the area of founding theory 
for the FA S method. 

A new technique, called multilevel projection methods (PML) has recently been introduced, [18] 
in an effort to provide a unifying, thematic approach to the design of a multilevel solver for a given 
problem. The main feature of PML methods is that the only basic choices that must be made 
concern the subspaces that will be used in relaxation and coarsening. All other components of the 
method, such as interlevel transfers, scaling, coarse-level problems, etc., are determined by 
projection between appropriate subspaces. In [18], several prototypical problems are developed to 
illustrate the principals involved. It now appears that the best hope of obtaining a strong founding 
theory for multilevel treatment of nonlinear problems may well be through careful and judicious 
application of PML , and our future research into solution methods for the problems we have 
discussed here will be aimed in that direction. 
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